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ARSTRACT  Environmental  geospatial  data  and  adult  and  larval  mosquito  collection  data  for  up  to 
106  sites  throughout  the  Republic  of  Korea  (ROK)  were  used  to  develop  ecological  niche  models 
(ENMs)  of  the  potential  geographic  distribution  for  eight  anopheline  species  known  to  occur  there. 
The  areas  predicted  suitable  for  the  Hyrcanus  Group  species  were  the  most  extensive  for  Anopheles 
sinensis  Wiedemann,  An.  kleini  Rueda,  An.  belenrae  Rueda,  and  An.  pullus  Yamada,  intermediate  for 
An.  sineroides  Yamada,  and  the  most  restricted  for  An.  lesteri  Raisas  and  Hu  and  the  non-Hyrcanus 
Group  species  An.  koreicus  Yamada  and  Watanabe  and  An.  lindesayi  Yamada.  The  relative  vectorial 
importance  of  these  species  is  unknown,  and  all,  except  An.  koreicus  and  An.  lindesayi,  are  predicted 
to  occur  widely  in  the  northwest  of  the  ROK  where  malaria  transmission  has  been  sporadic  since  its 
resurgence  in  1993.  Our  ENMs  suggest  that  it  is  unlikely  that  An.  koreicus  and  An.  lindesayi  are  vectors, 
but  we  do  not  document  consistent  geographic  differentiation  that  might  incriminate  any  of  the  other 
species  as  vectors,  because  all  species  are  predicted  to  occur  in  North  Korea,  we  also  cannot  reject 
the  hypothesis  that  malaria  infected  mosquitoes  from  North  Korea  may  have  been  the  cause  of  the 
resurgence  of  malaria  in  the  ROK.  Ecological  differentiation  of  the  eight  species  is  inferred  from 
collection  locations  and  34  environmental  layers  based  on  remote  sensing  and  global  climatic  averages. 
Interspecific  differences  were  noted,  and  characterizing  mosquito  habitats  by  ground-based  and 
remote  sensing  methods  is  proposed. 
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Our  understanding  of  malaria  transmission  risks  and 
the  types  of  vector  control  needed  for  any  given  lo¬ 
cation  relies  in  part  on  what  is  known,  or  assumed, 
about  the  geographical  distribution  of  each  of  the 
potential  vector  species.  The  availability  of  accurate 
distribution  maps  for  vector  species,  in  combination 
with  demographic  and  malaria  diagnosis  data,  could 
assist  in  determining  the  extent  of  the  “Mal-area” 
where  malaria  transmission  is  possible  (Foley  et  al. 
2008a). 

The  recent  availability  of  tools  for  ecological  niche 
modeling  (ENM)  has  made  it  possible  to  outline  the 
ecological  niche  and  potential  geography  of  a  variety 
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of  mosquito  vectors  (Levine  et  al.  2004,  Moffett  et  al. 
2007,  Foley  et  al.  2008b).  Ecological  niches  are  mod¬ 
eled  using  point-occurrence  data  as  they  relate  to 
digital  geospatial  data  layers  summarizing  environ¬ 
mental  variation.  Ruermann  et  al.  (2008)  found  that  a 
combination  of  remote  sensing  and  climatic  layers 
resulted  in  better  ENM  performance  than  using  these 
layers  separately.  We  were  interested  in  using  remote 
sensing  and  climatic  layers  to  model  the  distribution  of 
potential  malaria  vectors  in  the  Republic  of  Korea 
(South  Korea,  ROK).  Matching  mosquito  distribution 
to  foci  of  malaria  transmission  could  help  incriminate 
those  species  that  are  the  primary  malaria  vectors 
(Peterson  2006a). 

Eefore  the  1950s,  Plasmodium  vivax  malaria  was 
endemic  and  widespread  in  the  ROK  (Ree  2000)  and 
was  eventually  eradicated  in  1979  (WHO  1981).  In 
1993,  malaria  re-emerged  and  reached  a  peak  of  4,142 
cases  in  2000  before  falling  to  774  cases  in  2004  (Lee 
et  al.  2007).  Most  malaria  cases  seem  to  have  been 
contracted  near  the  demilitarized  zone  (DMZ)  that 
separates  North  Korea  (Democratic  Peoples  Republic 
of  Korea,  DPRK)  and  the  ROK  (Lee  et  al.  2007).  One 
hypothesis  to  explain  the  resurgence  of  malaria  in  the 
ROK  is  that  infected  mosquitoes  traveled  across  the 
DMZ  from  the  DPRK  (Chai  et  al.  1994,  Kho  et  al.  1999, 
Ree  2000) .  In  this  case,  the  vector  species  responsible 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2009  2.  REPORT  TYPE 

4.  TITLE  AND  SUBTITLE 

Geographic  Distribution  and  Ecology  of  Potential  Malaria  Vectors  in 
Republic  of  Korea 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Walter  Reed  Army  Institute  of  Research, Walter  Reed  Biosystematics 
Unit, 503  Robert  Grant  Avenue, Silver  Spring, MD ,20910 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

Environmental  geospatial  data  and  adult  and  larval  mosquito  collection  data  for  up  to  106  sites  throughout 
the  Republic  of  Korea  (ROK)  were  used  to  develop  ecological  niche  models  (ENMs)  of  the  potential 
geographic  distribution  for  eight  anopheline  species  known  to  occur  there.  The  areas  predicted  suitable  for 
the  Hyrcanus  Group  species  were  the  most  extensive  for  Anopheles  sinensis  Wiedemann,  An.  kleini  Rueda, 
An.  belenrae  Rueda,  and  An.  pullus  Yamada,  intermediate  for  An.  sineroides  Yamada,  and  the  most 
restricted  for  An.  lesteri  Baisas  and  Hu  and  the  non-Hyrcanus  Group  species  An.  koreicus  Yamada  and 
Watanabe  and  An.  lindesayi  Yamada.  The  relative  vectorial  importance  of  these  species  is  unknown,  and 
all,  except  An.  koreicus  and  An.  lindesayi,  are  predicted  to  occur  widely  in  the  northwest  of  the  ROK 
where  malaria  transmission  has  been  sporadic  since  its  resurgence  in  1993.  Our  ENMs  suggest  that  it  is 
unlikely  that  An.  koreicus  and  An.  lindesayi  are  vectors  but  we  do  not  document  consistent  geographic 
differentiation  that  might  incriminate  any  of  the  other  species  as  vectors.  Because  all  species  are  predicted 
to  occur  in  North  Korea,  we  also  cannot  reject  the  hypothesis  that  malaria  infected  mosquitoes  from  North 
Korea  may  have  been  the  cause  of  the  resurgence  of  malaria  in  the  ROK.  Ecological  differentiation  of  the 
eight  species  is  inferred  from  collection  locations  and  34  environmental  layers  based  on  remote  sensing  and 
global  climatic  averages.  Interspeci?  differences  were  noted,  and  characterizing  mosquito  habitats  by 
ground-based  and  remote  sensing  methods  is  proposed. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

13 

3.  DATES  COVERED 

00-00-2009  to  00-00-2009 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


May  2009 


Foley  et  al.:  Distribution  of  Malaria  Vectors  in  Korea 


681 


for  the  resurgence  of  malaria  would  likely  be  present 
on  both  sides  of  the  DMZ  and  be  able  to  traverse  this 
environment,  possibly  feeding  off  nonhuman  hosts.  A 
recent  competing  hypothesis  is  that  malaria  arose 
spontaneously  in  both  South  and  North  Korea  as  a 
result  of  a  combination  of  troop  build-up  along  the 
DMZ,  increased  exposure  to  local  mosquitoes,  and  the 
occasional  relapse  caused  by  hypnozoites  (Hulden 
and  Hulden  2008). 

The  anopheline  fauna  of  the  ROK  is  relatively  well 
resolved  taxonomically  ( Wilkerson  et  al.  2003,  Li  et  al. 
2005,  Rueda  2005)  and  includes  six  species  of  the 
Hyrcanus  Group:  Anopheles  sinensis  s.s.  Wiedemann; 
An.  pullus  M.  Yamada  ( =An.  yatsushiroensis ) ;  An.  les- 
teri  Raisas  and  Hu  (—An.  anthropophagus );  An.  siner¬ 
oides  S.  Yamada;  An.  kleini  Rueda;  An.  belenrae  Rueda; 
and  two  non-Hyrcanus  Group  species  An.  lindesayi 
japonicus  S.  Yamada  and  An.  koreicus  S.  Yamada  and 
Watanabe.  Not  all  of  these  species  can  be  identified  on 
morphology,  and  a  polymerase  chain  reaction  (PCR) 
technique  has  been  developed  for  species  identifica¬ 
tion  (Li  et  al.  2005).  Some  information  has  been  col¬ 
lected  on  the  distribution  and  larval  habitat  charac¬ 
teristics  of  molecularly  identified  anophelines  of  the 
ROK  (Rueda  et  al.  2006,  Kim  et  al.  2007,  Sanies  et  al. 
2008).  Previous  attempts  to  incriminate  the  malaria 
vector  species  are  flawed  because  of  an  imperfect 
understanding  of  the  anopheline  fauna  in  the  ROK. 
Historically,  An.  sinensis  was  considered  the  primary 
vector.  However,  the  discovery  of  additional  species 
and  results  from  field  and  laboratory  parasite  studies 
have  suggested  that  An.  kleini  and  An.  pullus  are  the 
primary  vectors  and  An.  sinensis  is  a  secondary  vector 
(Lee  et  al.  2007) .  However,  because  further  mosquito 
and  parasite  sampling  is  needed,  we  regard  here  all 
species  as  candidate  vectors. 

Our  aim  is  to  present  ENM-based  distributional 
predictions  for  the  total  anopheline  species  known 
from  the  ROK  rather  than  for  any  subset  of  species.  A 
multiple  species  approach  is  potentially  more  infor¬ 
mative  than  a  single  species  approach  but  has  not  been 
commonly  undertaken  until  recently  (Rarros  et  al. 
2007,  Moffett  et  al.  2007,  Rosa-Freitas  et  al.  2007) .  The 
resulting  models  of  potential  species  distributions  can 
be  compared  with  areas  of  suspected  malaria  trans¬ 
mission  to  better  understand  the  historical  and  current 
distribution  of  malaria  risks  in  the  ROK.  In  particular, 
we  are  interested  in  determining  whether  geograph¬ 
ical  differentiation  might  incriminate  any  one  species 
as  a  vector  or  the  species  responsible  for  malaria  re¬ 
surgence.  The  ecological  requirements  of  each  species 
were  studied  using  data  derived  from  remote  sensing 
and  global  climatic  averages  available  for  each  collec¬ 
tion  location  to  provide  further  insight  into  species 
distribution. 

Materials  and  Methods 

Species  Occurrence  Data.  The  occurrence  points 
comprised  152  unique  sites  for  An.  sinensis,  51  for  An. 
pullus,  54  for  An.  kleini,  49  for  An.  belenrae,  27  for  An. 
sineroides,  13  for  An.  koreicus,  16  for  An.  lesteri,  and  10 


for  An.  lindesayi.  These  were  for  larval  and  adult  col¬ 
lections  from  1998  to  2006,  some  of  which  have  been 
published  (Rueda  et  al.  2006,  Sames  et  al.  2008).  Spe¬ 
cies  were  identified  using  PCR  as  described  by  Li  et  al. 
(2005).  Clustering  of  sample  points  was  evident,  with 
most  occurring  in  the  northwest  of  the  ROK,  near  the 
DMZ.  To  minimize  the  effects  of  spatial  autocorrela¬ 
tion,  only  localities  separated  by  at  least  5  km  were 
used,  because  it  was  considered  that  these  sites  would 
exhibit  sufficient  potential  variation  to  be  considered 
spatially  independent.  Thus,  only  a  total  of  106  points 
were  used  in  the  final  analysis,  i.e.,  80  (An.  sinensis ) ,  33 
(An.  pullus),  32  (An.  kleini),  36  (An.  belenrae),  20  (An. 
sineroides) ,  9  (An.  koreicus) ,  12  (An.  lesteri),  and  9  (An. 
lindesayi).  For  the  purposes  of  model  evaluation  (see 
below),  separate  analyses  were  done  for  species  with 
low  numbers  of  sample  points  (<25)  and  those  with 
high  numbers  (£25).  All  occurrence  locality  data 
used  in  this  study  are  available  at  http://www. 
mosquitomap.org/  with  the  global  unique  identifi¬ 
ers  starting  with  MMap:KS. 

Environmental  Data.  We  obtained  GIS  raster  data 
layers  summarizing  elevation  (above  sea  level)  and  a 
selection  of  30  arc  second  resolution  climate  layers  for 
1950-1990  from  the  WorldClim  dataset  (Hijmans  et  al. 
2005,  http://www.worldclim.org/).  Climate  layers  in¬ 
cluded  monthly  temperature  (maximum,  minimum, 
mean),  monthly  precipitation,  and  19  "bioclimatic” 
layers  (representing  annual  trends,  seasonality,  and 
extreme  or  limiting  environmental  factors) .  These  lay¬ 
ers  were:  annual  mean  temperature,  mean  diurnal 
range  [mean  of  monthly  (max  temperature  —  min 
temperature)  =  P2],  isothermality  [(P2/P7)  X  100], 
temperature  seasonality  (SD  X  100),  maximum  tem¬ 
perature  of  the  warmest  month  (P5),  minimum  tem¬ 
perature  of  coldest  month  (P6),  temperature  annual 
range  [  (P5  —  P6)  =  P7],  mean  temperature  of  the 
wettest  quarter,  mean  temperature  of  the  driest  quar¬ 
ter,  mean  temperature  of  the  warmest  quarter,  mean 
temperature  of  the  coldest  quarter,  annual  precipita¬ 
tion,  precipitation  of  the  wettest  month,  precipitation 
of  the  driest  month,  precipitation  seasonality  (coeffi¬ 
cient  of  variation),  precipitation  of  the  wettest  quar¬ 
ter,  precipitation  of  the  driest  quarter,  precipitation  of 
the  warmest  quarter,  and  precipitation  of  the  coldest 
quarter.  Temperature  is  in  degrees  Celcius  (X10)  and 
precipitation  in  millimeters. 

We  also  obtained  five  layers  summarizing  topogra¬ 
phy  and  landform  (compound  topographic  index, 
slope,  aspect,  flow  direction,  and  flow  accumulation) 
from  the  U.S.  Geological  Survey’s  Hydro-IK  Elevation 
Derivative  Database  (http://edc.usgs.gov/products/ 
elevation /gtopo30/ hydro/,  native  resolution  1  by  1 
km).  The  compound  topographic  index  layer  summa¬ 
rizes  upward  curvature  of  land  surface  and  conse¬ 
quent  water  pooling,  which  is  important  for  mosquito 
development.  Aspect,  measured  in  degrees,  can  be 
thought  of  as  slope  direction  with  zero  at  north  and 
increasing  clockwise.  Flow  accumulation  at  this  res¬ 
olution  translates  directly  into  upstream  drainage  area 
in  square  kilometers. 
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We  included  data  layers  summarizing  the  monthly 
photosynthetic  capacity  or  “greenness  index”  termed 
the  normalized  difference  vegetation  index  (NDVI; 
0.0091°  resolution)  (Tucker  1979)  from  the  advanced 
very  high  resolution  radiometer  (AVHRR)  satellite 
(http://glcf.umiacs.umd.edu/ data/landcover/ )  for 
1992-1993.  Data  layers  summarizing  percent  tree 
cover  and  13  classes  of  land  use/land  cover  (0.0083° 
resolution)  were  obtained  from  the  Global  Land 
Cover  Facility  (http://glcf.umiacs.umd.edu/data/), 
derived  from  data  for  1992-1993  and  1981-1994,  re¬ 
spectively.  For  the  ecological  comparisons,  we  in¬ 
cluded  data  on  soil  taxonomy  suborders  (0.0333°  res¬ 
olution),  obtained  from  the  U.S.  Department  of 
Agriculture  National  Soils  Conservation  Service  (http:/  / 
soils.usda.gov/ use/ worldsoils/  mapindex/  order.html) . 
Udults  (Order:  Ultisols)  and  Orthents  (Order:  Enti- 
sols)  are  the  most  common  soil  types  in  the  study  area. 
According  to  Anon.  (1999),  Udults  are  more  or  less 
freely  drained,  relatively  humus-poor  Ultisols  that 
have  an  udic  (humid)  moisture  regimen.  Udults  occur 
in  humid  climates,  usually  receive  well  distributed 
rainfall,  are  often  forested  in  a  natural  state,  and  are 
often  cultivated.  Orthents  are  commonly  on  recent 
erosional  surfaces,  which  may  be  because  of  geological 
processes  or  induced  by  cultivation,  mining,  or  other 
factors.  Orthents  occur  in  any  climate  and  under  any 
vegetation  but  usually  do  not  occur  in  areas  that  have 
aquic  (saturated)  conditions,  as  occurs  with  a  high 
water  table. 

All  environmental  data  layers  were  resampled  with 
bilinear  interpolation  (Phillips  et  al.  2006)  to  a  pixel 
resolution  of 0.0083°  or  1  by  1  km.  Layers  were  clipped 
to  an  area  encompassing  the  area  of  interest  on  the 
Korean  peninsula  (Fig.  1;  33.0877-39.4875°  N, 
124.3686-129.885°  E),  for  a  total  of  508,416  pixels  per 
layer.  Recause  overfitting  is  more  likely  in  highly  di¬ 
mensional  environmental  spaces  (Peterson  et  al. 
2007a),  we  used  a  principal  components  analysis 
(PCA;  based  on  correlation  matrix;  in  Minitab  15.1.1.; 
Minitab,  State  College,  PA)  of  the  86  environmental 
layers  for  ENMs  to  create  15  new  axes  summarizing 
overall  variation  in  fewer,  independent  dimensions. 
Distribution  models  are  available  at  http://www. 
mosquitomap.org/ . 

ENMs.  We  used  the  genetic  algorithm  for  rule-set 
prediction  (Garp;  Stockwell  and  Noble  1992)  and  a 
maximum  entropy  approach  (Maxent;  Phillips  et  al. 
2006)  for  ENM  development.  These  algorithms  were 
selected  because  both  have  been  used  in  ENM  devel¬ 
opment  for  mosquitoes,  and  both  have  seen  intensive 
analysis  of  performance  (Peterson  et  al.  2007a,  Phillips 
2008). 

We  used  a  desktop  version  of  Garp  (Desktop  Garp 
1.1.3;  Scachetti-Pereira  2003).  Garp  uses  an  iterative 
process  of  rule  selection,  evaluation,  testing,  and  in¬ 
corporation  or  rejection.  A  rule  is  selected  and  applied 
to  the  training  data  (50-75%  of  points  input  into  the 
program).  The  genetic  algorithm  in  Garp  allows  rules 
to  “evolve”  to  maximize  predictive  accuracy.  The 
change  in  predictivity  between  iterations  is  used  to 
evaluate  whether  a  particular  rule  should  be  incorpo¬ 


rated  into  the  model,  and  the  algorithm  runs  until 
convergence  (set  at  0.01)  or  to  a  maximum  of  1,000 
iterations.  We  generated  up  to  1,000  models  for  each 
species  and  selected  the  10  “best”  models  based  on 
omission  and  commission  error  statistics  (Anderson  et 
al.  2003):  specifically,  we  used  an  absolute  omission 
threshold  of  10%  and  the  central  50%  of  the  commis¬ 
sion  index  distribution. 

Maxent  has  been  described  and  tested  in  detail 
(Phillips  et  al.  2006,  Phillips  and  Dudik  2008),  and  we 
used  version  3.2.1  (Phillips  2006).  This  version  has  a 
number  of  improvements,  including  a  default  logistic 
output  format,  which  can  be  interpreted  as  the  prob¬ 
ability  that  the  species  is  present,  conditional  on  the 
environmental  conditions  (Phillips  and  Dudik  2008). 
Maxent  is  based  on  the  idea  that  the  best  explanation 
for  unknown  phenomena  will  maximize  the  entropy  of 
the  probability  distribution,  subject  to  the  constraint 
of  the  environmental  conditions  observed  at  sites 
where  species  have  been  detected.  We  used  default 
parameters  for  Maxent  models,  except  that  a  random 
seed  was  used. 

Model  Evaluation.  For  Garp  output  of  the  four  spe¬ 
cies  with  25  or  more  sample  points,  extrinsic  test  data 
were  overlaid  onto  predicted  distributions  and  ob¬ 
served  correct  predictions  were  tallied.  The  observed 
degree  of  coincidence  was  tested  according  to  random 
expectations,  calculated  as  the  product  of  the  number 
of  testing  points  and  the  proportional  coverage  of  the 
study  area  predicted  to  be  present.  We  followed 
Anderson  et  al.  (2002)  and  Peterson  et  al.  (2007b)  in 
using  cumulative  binomial  probabilities  that  the  ob¬ 
served  success  could  be  the  result  of  random  associ¬ 
ation  between  predictions  and  test  points  as  a  test 
statistic.  For  Maxent  output  of  the  four  species  with  25 
or  more  sample  points,  model  validation  statistics  were 
based  on  10  replicate  random  partitions  of  the  local¬ 
ities  into  the  test  (25%)  and  training  (75%)  data.  For 
each  replicate,  we  calculated  the  number  of  test  lo¬ 
calities  omitted  from  the  prediction  and  applied  a 
binomial  test  to  check  statistical  significance  with  re¬ 
gard  to  the  proportion  of  the  study  area  that  was 
predicted  as  present. 

For  the  four  species  with  <25  sample  sites,  for  both 
Garp  and  Maxent  analyses,  we  used  the  jackknife  (or 
“leave-one-out”)  validation  approach  developed  by 
Pearson  et  al.  (2007),  using  the  program  made  avail¬ 
able  by  those  authors.  The  more  customary  receiver 
operating  characteristic  analysis  was  not  used  in  this 
study  because  of  recent  criticism  and  concerns  about 
bias  and  artifact  (Lobo  et  al.  2008,  Peterson  et  al. 
2008). 

Model  validation  and  interpretation  requires  setting 
a  decision  threshold  to  distinguish  areas  of  “presence” 
and  “absence.”  We  follow  the  approach  of  Pearson  et 
al.  (2007)  by  using  the  lowest  presence  threshold 
(LPT).  The  LPT  approach  is  relatively  conservative, 
identifying  pixels  that  are  at  least  as  suitable  as  the 
lowest  value  associated  with  a  species’  presence 
(Pearson  et  al.  2007). 

Consensus  Predictions.  Recause  some  studies 
(Peterson  et  al.  2007a)  have  suggested  that  Maxent 
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Fig.  1.  Predicted  lowest  presence  threshold  distribution  of  An.  belenrae.  An.  kleini.  An.  koreicus,  An.  lesteri,  An.  lindesayi, 
An.  pullus,  An.  sinensis,  and  An.  sineroides  in  Korea  using  the  genetic  algorithm  for  rule-set  prediction  (Garp).  Location  of 
known  occurrences  of  each  species  is  shown  by  black  dots. 


models  may  be  better  for  interpolation  and  Garp 
models  better  for  extrapolation  (i.e.,  transferabil¬ 
ity),  we  calculated  a  consensus  of  the  two  model 


outputs  for  analyses  trained  on  the  full  data  set  of 
sample  points  for  each  species.  Following  Foley  et 
al.  (2008b) ,  we  calculated  the  area  of  intersection  of 
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Garp  and  Maxent  models,  but  using  the  LPT  as  the 
decision  threshold. 

Variable  Importance  and  Response  Curves.  We 
used  the  program  PAST  1.27  (http://folk.uio.no/ 
ohammer/past/ )  to  undertake  a  nonmetric  multidi¬ 
mensional  scaling  analysis  of  the  Euclidian  distance  of 
means  of  the  first  15  principal  components  summa¬ 
rizing  the  86  environmental  data  layers  occurring  at 
Korean  collection  locations  for  the  eight  mosquito 
species.  This  analysis  is  appropriate  for  the  reduction 
and  interpretation  of  large  multivariate  ecological  data 
sets  and  was  done  to  assist  in  understanding  differ¬ 
ences  in  the  ecological  requirements  of  these  species. 
The  algorithm  in  PAST  1.27  attempts  to  place  the  data 
points  in  a  two-dimensional  coordinate  system  such 
that  the  ranked  differences  are  preserved. 

We  used  Maxent  to  investigate  response  curves  and 
variable  importance  using  33  environmental  variables 
(bioclimatic,  topography  and  landform,  NDVI,  tree 
cover,  land  cover,  and  elevation;  Table  1). 

Results 

The  results  from  the  PCA  of  environmental  data 
showed  that  the  first  principal  component  (PCI)  was 
most  closely  related  to  aspects  of  temperature, 
whereas  PC2  was  related  mainly  to  precipitation  and 
seasonality  and  PC3  to  precipitation,  land  cover,  and 
NDVI  (data  not  shown) .  Subsequent  PCs  emphasized 
NDVI  and  temperature  and  precipitation  extremes 
(PC4),  isothermality  (PC5),  land  cover  and  topogra¬ 
phy  (PC6  and  PC12),  precipitation  and  topography 
(PC7),  topography  (PC8  and  PC9),  precipitation 
(PC10  and  PC11),  precipitation  and  NDVI  (PC13), 
slope  and  NDVI  (PC14),  and  topography  (PC15).  In 
all,  the  first  15  PCs  captured  97.4%  of  the  overall 
variation  in  the  original  86  environmental  data  layers. 

Maxent:  Species  With  High  Numbers  of  Sample 
Points.  For  Maxent  predictions  of  An.  sinensis  using 
total  point  for  training,  the  LPT  score  was  0.0683.  For 
replicates  where  points  were  divided  into  test  and 
training  points.  An.  sinensis  omission  error  ranged 
from  0  to  0.1  (mean,  0.02;  SD,  0.035) ,  equivalent  to  0  -2 
test  from  20  test  localities  omitted,  with  8/ 10  partitions 
recording  minimum  training  presence  binomial  prob¬ 
ability  of  P  <  0.05.  For  Maxent  predictions  of  An. 
belenrae  using  total  point  for  training,  the  LPT  score 
was  0.1452.  For  replicates  where  points  were  divided 
into  test  and  training  points,  An.  belenrae  omission 
error  ranged  from  0  to  0.444  (mean,  0.178;  SD,  0.167), 
equivalent  to  0  -  4  test  from  nine  test  localities  omitted, 
with  6/10  partitions  recording  minimum  training  pres¬ 
ence  binomial  probability  of  P  <  0.05.  For  Maxent 
predictions  of  An.  pullus  using  total  point  for  training, 
the  LPT  score  was  0.1327.  For  replicates  where  points 
were  divided  into  test  and  training  points,  An.  pullus 
omission  error  ranged  from  0  to  0.625  (mean,  0.225; 
SD,  0.194),  equivalent  to  0-5  test  from  eight  test  lo¬ 
calities  omitted,  with  5/10  partitions  recording  mini¬ 
mum  training  presence  binomial  probability  of  P  < 
0.05.  For  Maxent  predictions  of  An.  kleini  using  total 
point  for  training,  the  LPT  score  was  0.1117.  For  rep¬ 


licates  where  points  were  divided  into  test  and  train¬ 
ing  points.  An.  kleini  omission  error  ranged  from  0  to 
0.429  (mean,  0.214;  SD,  0.168),  equivalent  to  0-4  test 
from  seven  test  localities  omitted,  with  4/ 10  partitions 
recording  minimum  training  presence  binomial  prob¬ 
ability  of  P  <  0.05.  Thus,  Maxent  distribution  models 
for  species  with  larger  sample  sizes  (25  or  more) 
indicated  that  An.  sinensis  had  the  best  support,  ac¬ 
cording  to  random  expectations  of  omission  error  and 
An.  kleini  the  least  support. 

Maxent:  Species  With  Low  Numbers  of  Sample 
Points.  For  Maxent  predictions  of  An.  sineroides  using 
total  point  for  training,  the  LPT  score  was  0.2668.  After 
Jackknife  validation,  An.  sineroides  omission  error  oc¬ 
curred  in  6  of  20  replicates,  with  average  minimum 
training  presence  area  =  0.412  and  P  =  0.005696.  For 
Maxent  predictions  of  An.  lindesayi  using  total  point 
for  training,  the  LPT  score  was  0.3982.  After  Jackknife 
validation,  An.  lindesayi  omission  error  occurred  in 
four  of  nine  replicates,  with  average  minimum  training 
presence  area  =  0.215  and  P  =  0.011035.  For  Maxent 
predictions  of  An.  lesteri  using  total  point  for  training, 
the  LPT  score  was  0.2412.  After  Jackknife  validation, 
An.  lesteri  omission  error  occurred  in  2  of  12  replicates, 
with  average  minimum  training  presence  area  =  0.034 
and  P  =  0.833333.  For  Maxent  predictions  of  An.  ko- 
reicus  using  total  point  for  training,  the  LPT  score  was 
0.3992.  After  Jackknife  validation,  An.  koreicus  omis¬ 
sion  error  occurred  in  2  of  10  replicates,  with  average 
minimum  training  presence  area  =  0.163  and  P  = 
0.000015.  Thus,  Maxent  distribution  models  for  species 
with  smaller  sample  sizes  (<25)  indicated  that  all 
species,  with  the  exception  of  An.  lesteri,  gave  better 
omission  error  than  random  expectations. 

Garp:  Species  With  High  Numbers  of  Sample 
Points.  For  Garp  predictions  of  An.  sinensis  using  total 
point,  the  LPT  score  was  threshold  5  (T5).  In  other 
words,  a  collection  point  was  not  predicted  as  present 
unless  at  least  5  of  10  of  the  best  models  agreed.  All 
thresholds  were  better  than  random  expectations  (P  < 
0.05) .  For  Garp  predictions  of  An.  kleini,  the  LPT  score 
was  T7.  All  thresholds  were  better  than  random  ex¬ 
pectations  (P  <  0.05).  For  Garp  predictions  of  An. 
belenrae,  the  LPT  score  was  T8.  All  thresholds  were 
better  than  random  expectations  (P  <  0.05) .  For  Garp 
predictions  of  An.  pullus,  the  LPT  score  was  T8.  All 
thresholds  were  better  than  random  expectations  (P  < 
0.05).  Thus,  Garp  distribution  models  for  species  with 
larger  sample  sizes  (25  or  more)  indicated  that  all 
species  gave  better  omission  error  than  random  ex¬ 
pectations. 

Garp:  Species  With  Low  Numbers  of  Sample  Points. 
For  Garp  predictions  of  An.  sineroides,  the  LPT  score 
T10.  After  Jackknife  validation,  An.  sineroides  omission 
error  occurred  in  17  of  20  replicates,  with  average 
minimum  training  presence  area  =  0.143799  and  P  = 
0.427249.  For  Garp  predictions  of  An.  lindesayi,  the 
LPT  score  T10.  After  Jackknife  validation,  An.  linde¬ 
sayi  omission  error  occurred  in  10  of  10  replicates,  with 
average  minimum  training  presence  area  =  0.0002  and 
P  =  1.0.  For  Garp  predictions  of  An.  lesteri,  LPT  was 
T10.  After  Jackknife  validation,  An.  lesteri  omission 
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error  occurred  in  7  of  11  replicates,  with  average 
minimum  training  presence  area  =  0.016503  and  P  — 
0.000007.  For  Garp  predictions  of  An.  koreicus,  the  LPT 
score  was  T10.  After  Jackknife  validation,  An.  koreicus 
omission  error  occurred  in  8  of  10  replicates,  with 
average  minimum  training  presence  area  =  0.002  and 
P  =  0.000009.  Thus,  Garp  distribution  models  for  spe¬ 
cies  with  smaller  sample  sizes  (<25)  indicated  that 
only  An.  lesteri  and  An.  koreicus  gave  better  omission 
error  than  random  expectations. 

Species  models  based  on  all  data  points  (Figs.  1  and 
2)  identified  distribution  patterns  that  differed  be¬ 
tween  species  and  algorithms.  Garp  predictions 
tended  to  be  more  extensive  than  Maxent  predictions, 
as  has  been  found  in  other  studies  (Peterson  et  al. 
2007a) .  The  predicted  potential  distribution  of  species 
more  commonly  collected  (An.  belenrae,An.  kleini,An. 
pullus,  and  An.  sinensis)  showed  the  most  extensive 
distribution.  The  predicted  potential  distribution  of 
species  least  commonly  collected  (An.  koreicus.  An. 
lesteri,  An.  lindesayi,  and  An.  sineroides )  showed  the 
least  extensive  distribution.  Garp  and  Maxent  models 
differed  most  noticeably  for  An.  koreicus  and  An.  linde¬ 
sayi,  which  also  had  the  lowest  numbers  of  sample 
points.  The  Maxent  model  for  An.  koreicus  identified 
three  distinct  areas  of  occurrence  in  the  ROK  com¬ 
pared  with  one  diffuse  area  in  the  Garp  model.  The 
Maxent  model  for  An.  lindesayi  predicts  that  this  spe¬ 
cies  occurs  more  to  the  south  in  the  ROK  than  the 
Garp  model  (which  did  not  significantly  differ  from 
random  expectation).  The  Maxent  model  for  An.  sin¬ 
eroides  predicts  that  this  species  occurs  along  the  east 
coast  and  has  a  different  southern  distribution  in  the 
ROK  compared  with  the  Garp  prediction  (which  did 
not  significantly  differ  from  random  expectation) .  The 
intersection  of  the  models  using  the  LPT  threshold  is 
shown  in  Fig.  3  for  an  area  north  of  Seoul,  where 
resurgent  malaria  transmission  has  occurred.  All  spe¬ 
cies  were  predicted  to  occur  in  the  northwest  of  the 
ROK  and  in  North  Korea. 

Variable  Importance  and  Response  Curves.  Output 
distribution  models  using  bioclimatic,  topography  and 
landform,  NDVI,  tree  cover,  land  cover,  and  elevation 
variables  (see  Table  1 )  were  similar  to  but  less  general 
than  models  using  PCA-derived  variables.  Analysis  of 
the  jackknife  of  regularized  training  gain  showed  that 
the  NDVI  for  June  gave  the  most  useful  unique  in¬ 
formation  (i.e.,  decreased  the  gain  the  most  when  it 
was  omitted)  for  most  species  and  the  precipitation  of 
the  wettest  month  resulted  in  the  highest  gains  when 
used  in  isolation  (i.e.,  the  most  useful  information  by 
itself)  for  four  of  the  eight  species. 

Figure  4  shows  a  minimal  spanning  tree  of  the  non¬ 
metric  multidimensional  scaling  analysis  of  the  means 
of  the  first  15  principal  components  for  collection 
locations  of  the  eight  mosquito  species.  Average  en¬ 
vironmental  conditions  for  the  non-Hyrcanus  Group 
species  An.  koreicus  and  An.  lindesayi  were  different  to 
Hyrcanus  Group  species.  However,  this  analysis 
showed  that  the  environment  of  An.  lesteri  was  the 
most  dissimilar  of  the  Hyrcanus  Group  species. 


Summary  data  for  the  eight  mosquito  species  for 
1-km2  grids  that  coincided  with  collection  locations  is 
given  in  Table  1  for  the  above  33  environmental  layers 
plus  soils.  Because  of  the  potential  for  bias  because  of 
some  species  (e.g.,  An.  sinensis)  being  more  exten¬ 
sively  distributed  in  the  warmer  southern  regions  of 
the  ROK  than  other  species,  comparison  of  environ¬ 
mental  conditions  at  locations  in  Table  1  were  limited 
to  those  occurring  north  of  latitude  36.5°  N  (all  species 
were  recorded  from  this  region) .  Based  on  these  data, 
the  following  environmental  conditions  were  charac¬ 
teristic  of  the  Hyrcanus  Group  species  (compared 
with  non-Hyrcanus  Group  species):  higher  annual 
mean  temperature,  lower  mean  diurnal  and  annual 
temperature  range,  lower  annual  precipitation  (par¬ 
ticularly  low  precipitation  in  the  warmest  and  wettest 
quarter),  lower  precipitation  seasonality,  lower  ele¬ 
vation,  lower  slope,  southerly  rather  than  south  west¬ 
erly  aspect,  lower  flow  accumulation  (i.e.,  less  up¬ 
stream  drainage),  lower  tree  cover  (except  for  An. 
sineroides),  occurring  in  cropland  and  wooded  grass¬ 
land  rather  than  deciduous  broadleaf  forest,  and  on 
Udults  soil  suborder  types  (except  for  An.  sineroides) 
rather  than  on  Orthents. 

The  following  lists  relative  tendencies  of  environ¬ 
mental  conditions  for  Korean  species  within  the  Hyr¬ 
canus  Group:  An.  sinensis,  higher  NDVI  in  most 
months  and  not  as  cold  in  the  coldest  part  of  the  year; 
An.  belenrae,  high  flow  accumulation  (i.e.,  greater  up¬ 
stream  drainage)  and  predominantly  cropland;  An. 
kleini,  southwesterly  aspect;  An.  pullus,  coldest  in  the 
coldest  part  of  the  year,  high  topoindex  (wetness), 
high  flow  accumulation  (i.e.,  greater  upstream  drain¬ 
age)  ,  low  slope,  and  bimodal  aspect;  An.  sineroides,  low 
temperature  seasonality  and  range  (not  as  cold  in  the 
coldest  part  of  the  year  and  not  as  warm  in  the  warmest 
part  of  the  year),  low  precipitation  seasonality  (not  as 
wet  in  the  wettest  part  of  the  year  and  not  as  dry  in  the 
driest  part  of  the  year) ,  low  topoindex  (wetness) ,  high 
slope,  low  flow  accumulation  (i.e.,  less  upstream 
drainage),  high  elevation,  high  tree  cover,  and  pre¬ 
dominantly  Orthents  soil  type;  An.  lesteri,  high  pre¬ 
cipitation  and  temperature  seasonality  (drier  and  wet¬ 
ter,  and  warmer),  low  mean  diurnal  temperature 
range,  low  isothermality,  low  NDVI  for  all  months,  low 
slope,  southwesterly  aspect,  low  tree  cover  (typically 
cropland) ,  and  low  elevation. 

A  one-way  analysis  of  variance  (ANOVA)  of  data 
for  each  environmental  layer  was  performed  using 
Minitab  for  the  Hyrcanus  Group  species.  This  analysis 
indicated  a  significant  difference  in  means  for  isother¬ 
mality  (F  =  3.60,  df  =  5,  P  =  0.004)  and  mean  tem¬ 
perature  of  the  wettest  quarter  (F  =  2.48,  df  =  5,  P  = 
0.034) ,  with  a  nearly  significant  result  for  mean  tem¬ 
perature  of  the  warmest  quarter  (F  =  2.26,  df  =  5,  P  = 
0.051 ).  In  all  cases  An.  lesteri  was  significantly  different 
from  the  other  species  (P  =  0.05) ,  based  on  the  Fisher 
least  significant  difference  (LSD)  method.  An.  lesteri 
has  low  scores  for  mean  diurnal  temperature  range  but 
high  scores  for  temperature  annual  range,  resulting  in 
low  scores  for  isothermality.  This  suggests  that  this 
species  prefers  warmer  and  more  stable  temperatures 
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Fig.  2.  Predicted  lowest  presence  threshold  distribution  of  An.  belenrae,  An.  kleini.  An.  koreicus.  An.  lesteri ,  An.  lindesayi. 
An.  pullus,  An.  sinensis,  and  An.  sineroides  in  Korea  using  a  maximum  entropy  (Maxent)  approach.  Location  of  known 
occurrences  of  each  species  is  shown  by  black  dots. 


during  the  summer  than  other  species  within  the  Hyr-  2.31,  df  =5 ,P  —  0.047)  for  PC5,  which  mainly  reflects 
canus  Group.  Only  one  significant  differences  for  the  the  differences  between  An.  lesteri  and  the  other  Hyr- 
15  PC  layers  was  recorded  (one-way  ANOVA,  F  =  canus  Group  species  in  terms  of  isothermality. 
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Fig.  3.  Consensus  by  intersection  of  predicted  lowest  presence  threshold  distribution  models  of  Garp  and  Maxent 
distributions  of  An.  belenrae.  An.  kleini.  An.  koreicus.  An.  lesteri,  An.  lindesayi.  An.  pullus.  An.  sinensis,  and  An.  sineroides.  The 
area  shown  encompasses  the  zone  of  recent  malaria  transmission  in  Korea.  Pixels  where  mosquito  species  are  predicted 
present  are  colored  gray  except  those  (in  black)  that  fall  within  areas  (outlined)  where  civilian  cases  of  vivax  malaria  were 
recorded  by  Kho  et  al.  (1999)  for  1994-1997.  Thus,  the  larger  the  amount  of  black  pixels  the  greater  the  evidence  for  that 
species  being  a  vector. 


Fig.  4.  Minimal  spanning  tree  of  a  nonmetric  multidi¬ 
mensional  scaling  analysis  of  the  means  of  the  first  15  prin¬ 
cipal  components  summarizing  86  environmental  data  layers 
occurring  at  Korean  collection  locations  of  An.  belenrae,  An. 
kleini,  An.  koreicus.  An.  lesteri.  An,  lindesayi,  An.  pullus,  An. 
sinensis,  and  An.  sineroides. 


Discussion 

The  species  identity  of  malaria  vectors  in  the  ROK 
is  unknown,  so  we  developed  ENMs  for  all  anopheline 
species  within  the  ROK  to  shed  light  on  their  distri¬ 
butions,  ecology,  and  possible  relation  to  malaria 
risk  in  that  country.  Prediction  successes  of  distri¬ 
bution  models  were  better  than  random  except  for 
Garp  models  for  An.  sineroides  and  An.  lindesayi  and 
Maxent  models  for  An.  lesteri.  An.  belenrae.  An.  pullus, 
and  An.  sinensis  have  been  variously  recorded  at  sites 
in  the  DPRK  by  Rueda  and  Gao  (2008)  (georeferenc- 
ing  error  noted  for  the  Pyeongseong  City  collection) . 
The  Kaesong  City  sites  were  positive  for  these  three 
species  in  our  predictions.  Although  our  Garp  predic¬ 
tion  for  An.  lindesayi  agrees  with  many  of  the  sites 
listed  in  Sanies  et  al.  (2008),  it  did  not  predict  any 
occurrence  on  Jeju  Island.  Lee  (1994)  reported  this 
species  on  Jeju  Island  (Korea),  but  this  location  may 
be  marginal,  because  subsequent  attempts  to  collect 
this  species  there  failed  (Kim  et  al.  2005,  Sames  et  al. 
2008).  The  spatial  resolution  of  this  analysis  (1  km2) 
may  may  be  too  coarse  for  smaller  habitats,  and  the 
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Aspect  (degrees  from  North) 


Fig.  5.  Frequency  distribution  of  aspect  (degrees  clock¬ 
wise  from  north)  for  Korean  collection  locations  for  three 
species  (An.  sinensis.  An.  pullus,  and  An.  lindesayi). 

number  of  locations  for  this  species  available  for  mod¬ 
eling  was  low.  Inclusion  of  additional  records,  espe¬ 
cially  at  lower  latitudes,  should  improve  predictions. 

One  hypothesis  to  explain  the  resurgence  of  malaria 
in  the  ROK  is  that  infected  mosquitoes  traveled  across 
the  DMZ  from  the  DPRK  (Chai  et  al.  1994,  Kho  et  al. 
1999,  Ree  2000).  Our  results  indicate  that  all  eight 
mosquito  species  considered  here  occur  on  both  sides 
of  the  DMZ  and  are  candidates  for  a  north  to  south 
border  crossing.  Thus,  we  did  not  document  geo¬ 
graphic  differentiation  that  could  reject  the  hypoth¬ 
esis  that  malaria-infected  mosquitoes  from  North  Ko¬ 
rea  may  have  been  the  cause  of  the  introduction  and 
resurgence  of  malaria  in  the  ROK. 

Did  we  document  geographic  differentiation  that 
might  incriminate  any  of  the  other  species  as  a  current 
or  past  vector?  Historically,  malaria  was  widespread  in 
the  ROK;  for  instance,  a  1965  map  from  the  National 
Malaria  Eradication  Service  in  Chai  (1999)  showed 
malaria  widespread  but  with  focal  areas  of  greater 
incidence  in  northern  Gyeongsangbuk  and  the  north 
and  southeast  of  Gyeonggi  provinces.  All  species  ex¬ 
cept  An.  lesteri,  An.  koreicus,  and  An.  lindesayi  are 
predicted  to  be  extensive  in  both  of  these  areas.  Mur¬ 
doch  and  Lueders  (1953)  reported  a  0.83%  malaria 
parasite  infection  rate  for  An.  sinensis  s.l.  collected  in 
the  Wonju  area  in  August  1952.  At  that  time,  morpho¬ 
logically  identified  An.  sinensis  could  include  other 
members  of  the  Hyrcanus  Group,  but  not  An.  siner- 
oides,An.  lindesayi,  and  An.  koreicus.  We  found  that  all 
species,  except  An.  koreicus.  An.  lesteri,  and  An.  linde¬ 
sayi,  are  predicted  to  occur  extensively  in  that  region. 
Regarding  more  recent  malaria  transmission.  Fig.  3 
shows  areas  where  civilian  cases  of  vivax  malaria  were 
recorded  in  1994-1997  (from  Kho  et  al.  1999).  All 
species  are  predicted  to  be  present  where  transmis¬ 
sion  is  thought  to  have  occurred,  but  the  least  co¬ 
occurrence  of  species  and  malaria  was  recorded  for 
An.  koreicus  and  An.  lindesayi,  with  An.  sineroides  dis¬ 
tribution  failing  to  co-occur  with  cases  in  the  south. 
These  species  are  also  infrequently  collected  in  adult 


and  larval  surveillance  (Kim  et  al.  2007).  Thus,  distri¬ 
bution  modeling  suggests  that  An.  koreicus  and  An. 
lindesayi,  at  least,  are  unlikely  to  have  played  a  sig¬ 
nificant  role  in  past  and  present  malaria  transmission 
in  the  ROK.  However,  it  is  possible  that  more  than  one 
anopheline  species  has  been  involved  in  malaria  trans¬ 
mission  in  the  ROK. 

Refore  eradication  in  the  1970s,  malaria  transmis¬ 
sion  occurred  widely  throughout  the  ROK,  possibly 
because  of  local  conditions,  such  as  agricultural  prac¬ 
tices,  favoring  vector  abundance,  and  survival,  or  be¬ 
cause  of  the  extent  of  human/ mosquito  contact  and 
lack  of  rapid  medical  treatment  that  favors  transmis¬ 
sion.  The  fact  that  resurgent  malaria  has  not  returned 
to  preeradication  limits  suggests  that  conditions  in  the 
ROK  have  changed.  For  example,  Korean  forests  were 
badly  degraded  through  the  first  half  of  the  20th  cen¬ 
tury,  but  between  1967  and  1995,  forested  areas  in¬ 
creased  to  6.3  million  ha  as  a  result  of  large-scale 
reforestation  (Food  and  Agriculture  Organization  of 
the  United  Nations  Forestry  county  profiles,  http:/ / 
www.fao.org/forestry/23809/en/kor/).  The  seasonal 
incidence  of  malaria  in  1913  versus  1993-1996  re¬ 
ported  in  Chai  (1999)  showed  a  shift,  from  June  and 
August  peaks  to  only  August,  which  could  indicate 
either  change  in  the  incubation  period  of  the  parasite, 
the  involvement  of  different  vector  species,  or  a 
change  in  the  phenology  of  the  vector  species.  Recent 
evidence  suggests  that  =60%  of  the  malaria  cases  are 
latent  with  the  onset  of  symptoms  6-18  mo  after  in¬ 
fection  (T.  A.  Klein,  personal  communication). 

Considering  ecological  requirements,  Hyrcanus 
Group  species  seemed  to  differ  from  non-Hyrcanus 
Group  species  (Fig.  4),  and  An.  lesteri  was  the  most 
divergent  among  the  former.  The  frequency  of  occur¬ 
rence  of  An.  sinensis  increases  to  a  maximum  for  lo¬ 
cations  with  a  south-facing  aspect,  whereas  aspect  is 
predominately  south  southeast  and  south  southwest 
for  An.  pullus  and  in  the  southwest  and  northwest 
quadrants  for  An.  lindesayi  (Fig.  5).  No  trend  with 
time  of  year  could  be  discerned  (data  not  shown). 
Could  this  difference  be  a  function  of  species  using 
different  land  cover  categories,  which  in  turn  occur 
more  frequently  at  different  aspects?  We  graphed  the 
aspect  of  all  grid  cells  (0.0083°  resolution)  for  Korea 
according  to  land  cover  category  (Fig.  6) .  Forest  areas 
(categories  1-5)  were  more  frequently  northeast  to 
southeast  facing,  whereas  crops  and  grassland  areas 
(7-11)  were  more  frequently  west  to  southwest  fac¬ 
ing.  The  forest  finding  agrees  with  Johnson  et  al. 
(2002),  who  showed  that  site  quality  for  oaks  in  West 
Virginia  (a  similar  latitude  to  the  ROK)  was  maximal 
at  81°  azimuth  and  least  favorable  at  261°;  increasing 
slope  reduced  site  quality  caused  by  associated  in¬ 
creases  in  solar  radiation.  The  cool  and  moist  micro¬ 
environment  in  the  northeast  quadrant  resulted  in 
thicker  leaf  litter  layers,  which  have  various  effects  on 
soil  development  processes  (Johnson  et  al.  2002). 
However,  despite  An.  lindesayi  and  An.  koreicus  being 
most  frequent  in  forest  areas,  they  are  more  likely  to 
be  found  in  areas  facing  southwest,  and  despite  An. 
sinensis  being  most  frequently  found  in  crops  and 
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Fig.  6.  Frequency  distribution  of  aspect  (degrees  clockwise  from  north)  for  grid  cells  (1  km2)  in  the  study  area  of  Korea, 
according  to  landscape  category  (1  =  evergreen  needleleaf  forest,  4  =  deciduous  broadleaf  forest,  5  =  mixed  forest,  6  = 
woodland,  7  =  wooded  grassland,  8  =  closed  shrubland,  9  =  open  shrubland,  10  =  grassland,  11  =  cropland,  12  =  bare  ground, 
and  13  =  urban  and  built) .  Note  that  forest  areas  (categories  1-5)  are  more  frequently  east  to  southeast  facing,  whereas  crops 
and  grassland  areas  (7-11)  are  more  frequently  west  to  southwest  facing.  Categories  2-3  are  not  present. 


grassland  areas,  it  is  most  often  found  in  land  facing 
south.  The  interaction  of  environmental  layers  may  be 
important  for  mosquito  ecology.  For  example,  slope 
and  aspect  would  be  expected  to  interact  with  tree 
cover  that  may  relate  to  a  preference  for  shaded  hab¬ 
itats.  As  mentioned  above,  land  cover  in  the  ROK  is 
dynamic.  Therefore,  landscape  and  tree  cover  esti¬ 
mates  from  remote  sensing  data  drawn  from  the  1980  - 
1990s  may  not  apply  to  mosquito  collections  con¬ 
ducted  after  2000. 

According  to  Kim  et  al.  (2007),  An.  kleini  numbers 
peak  during  late  spring  and  early  summer,  whereas  An. 
sinensis  peaks  during  late  summer  and  fall.  The  more 
southerly  extension  of  An.  sinensis  into  tropical  regions 
of  Indochina  supports  the  idea  that  this  species  is 
adapted  to  higher  temperatures.  Much  of  the  differ¬ 
ence  in  average  temperatures  for  An.  sinensis  and  An. 
kleini  is  caused  by  temperatures  during  the  colder 
months.  Overwintering  adaptations  such  as  diapause 
may  not  be  well  developed  in  An.  sinensis.  In  contrast, 
the  non-Hyrcanus  Group  species  were  found  in  loca¬ 
tions  with  the  lowest  (below  freezing)  winter  tem¬ 
peratures  and  may  have  well-developed  overwinter¬ 
ing  adaptations.  An.  lindesayi  have  been  found 
overwintering  in  the  ROK  as  larvae,  even  in  icy  water 
(H.-C.  Kim,  personal  communication). 

Kim  et  al.  (2007)  found  larvae  of  An.  sinensis  mainly 
from  water  seepage  springs  and  marshes,  An.  pullus 
from  rice  fields  and  stream  margins,  An.  belenrae  from 
pond  and  lake  habitats,  and  An.  lindesayi  from  pools  of 
fast  running  streams.  Some  concordance  between  re¬ 
sults  from  larval  habitat  studies  and  remote  sensing 
data  can  be  seen.  For  example,  we  found  that  An. 


pullus  occurred  mainly  in  croplands  with  low  slope, 
high  flow  accumulation,  and  high  wetness  (topoin- 
dex) ,  as  you  would  expect  of  rice  field  habitats,  and  An. 
lindesayi  occurred  in  areas  of  high  flow  accumulation 
and  high  slope,  as  you  would  expect  of  fast  running 
streams.  Combining  data  on  larval  microenvironment 
and  high-resolution  remotely  sensed  data  may  yield 
powerful  insights  into  ecological  differentiation. 

Maxent  response  variables  identified  NDVI  and  cli¬ 
mate  data  as  important  for  prediction  models,  thereby 
supporting  the  contention  of  Ruermann  et  al.  (2008) 
that  a  combination  of  remote  sensing  and  climatic 
layers  is  more  important  for  ENMs  than  either  is  in 
isolation.  Our  distribution  models  are  based  largely  on 
yearly  climate  averages  (e.g.,  WorldClim),  and  as 
such,  portray  general  potential  distributional  tenden¬ 
cies.  In  addition,  the  realism  of  ENMs  will  depend  in 
part  on  the  biological  relevance  for  mosquito  biology 
of  environmental  layers  such  as  NDVI  and  aspect. 
Absence  data  are  less  informative  for  assessing  model 
accuracy  because  ENM  shows  potential  distributions 
rather  than  actual  distributions;  other  factors  such  as 
historical,  physical,  climatological,  and  biotic  con¬ 
straints  may  limit  distribution  (Soberon  and  Peterson 
2005,  Peterson  2006b).  Improvements  in  model  accu¬ 
racy  could  be  made  by  using  more  widely  distributed 
sample  points,  absence  as  well  as  presence  informa¬ 
tion,  analyses  by  month  rather  than  year,  and  by  using 
more  up-to-date,  higher  resolution  environmental  lay¬ 
ers.  For  example,  more  detailed  information  on  rice 
irrigation,  tree  cover,  and  landscape  classes  may  im¬ 
prove  predictions. 
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Although  these  models  propose  to  identify  the  po¬ 
tential  limits  of  species  distributions,  they  do  not  de¬ 
fine  population  densities  and  survival,  two  critical 
components  of  vectorial  capacity.  For  example.  An. 
hleini  and  An.  pullns,  although  relatively  widely  dis¬ 
tributed,  in  agreement  with  model  predictions,  ac¬ 
counted  for  >50%  of  all  anophelines  collected  at  mil¬ 
itary  installations  near  the  DMZ.  South  of  Seoul,  these 
two  species  account  for  <5%  of  the  total  population  of 
mosquitoes  in  both  larval  and  adult  collections.  The 
distribution  of  malaria  transmission  in  the  ROK  has 
changed  in  the  last  century,  but  it  is  unknown  what 
role  changes  in  species  distribution,  if  any,  have  played 
on  changes  in  malaria  distribution.  Weather  anomalies 
and  disasters,  including  war,  could  affect  the  distri¬ 
bution  and  abundance  of  mosquitoes  because  of,  for 
instance,  changes  in  landscape  and/or  farming  prac¬ 
tices.  Thus,  mosquito  species  distribution  models  are 
only  the  first  step  in  the  process  of  modeling  the 
entomological  risk  factors  involved  in  malaria  trans¬ 
mission. 
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